Variational wavefunction study of the triangular lattice supersolid 
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We present a variational wavefunction which explains the behaviour of the supersolid state formed 
by hard-core bosons on the triangular lattice. The wavefunction is a linear superposition of only and 
all configurations minimising the repulsion between the bosons (which it thus implements as a hard 
constraint). Its properties can be evaluated exactly — in particular, the variational minimisation of 
the energy yields (i) the surprising and initially controversial spontaneous density deviation from 
half-filling (ii) a quantitatively accurate estimate of the corresponding density wave (solid) order 
parameter. 
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There has been sustained interest in using ultracold 
atoms confined in optical lattice potentials to realize 
strongly-correlated systems of interest to condensed mat- 
ter physics The recently discussed possibility that 
Helium has a supersolid phase Q leads, in this context, 
to a natural question: Can the lattice analog of this, 
namely a superfluid phase that simultaneously breaks lat- 
tice translation symmetry, be seen in such optical lattice 
experiments? 

Although other examples are known [2,], perhaps the 
best candidate for such a lattice supersolid persisting over 
a wide range of parameters (such as chemical potential 
and interaction energy) is that observed numerically in 
several recent Quantum Monte-Carlo (QMC) studies of 
a two-dimensional system of strongly interacting bosons 
in a triangular lattice potential .4j^], and confirmed in 
subsequent follow-up work Q . Following earlier work 0] , 
these QMC studies considered the model Hamiltonian: 



H = -tY,{b% + b]h) + V 



l/2)(n,-l/2)(l) 



where h\{hi) is the boson creation (annihilation) opera- 
tor, t represents the strength of the boson hopping am- 
plitude, V is the strength of the nearest neighbour re- 
pulsion, and the bosons are restricted to be in the hard- 
core limit {rii = 0,1) by a strong onsite repulsive term 
not written down explicitly (here, we have only displayed 
the Hamiltonian for the value of the chemical potential /i 
at which the system has particle- hole symmetry). Strik- 
ingly, in this particle-hole symmetric, hard-core case, an 
extended supersolid state was observed numerically for 
all V/t > 8.9, and seen to possess both a non-zero su- 
perfluid stiffness, and density wave ('solid') order at the 
three sublattice (-y/S x \/3 ordering) wavevector Q (the 
state is stable to changes in /i as well). 

There are two other genuinely surprisingly features of 
this phase. The first [3] concerns the nature of the solid 
order. Plausible mean-field theory arguments [H] predict 



that the density wave order in the supersolid at /i = 
should involve a '(H — 0)' type three sub-lattice pattern 
of density with pa = 1/2, pi, = 1/2 -\- 5, Pc — 1/2 — S, 

while in reality the system prefers a different '(H )' 

pattern with the same ordering wavevector: pa = 1/2 + 
h, Pb = Pc = 1/2 - So, with (5i,2 > (Fig[l|). General 
symmetry arguments 0, 0] predict that the total density 
of the system should also exhibit a spontaneous deviation 

from half-filling in a '(-I )' state. The second surprise 

is that although there is such a deviation, it is extremely 
small in 'natural' units (i.e compared to the strength of 
the density wave order itself), and therefore extremely 
difficult to see numerically 0, S] ■ 

In order to understand this large V/t supersolid phase 
in the hardcore limit, it is necessary to take into account 
this strong repulsion as a hard constraint. Here we con- 
sider an analytically tractable variational wavefunction 
that takes into account this hard constraint from the 
start, and uses a variational parameter to tunc the am- 
plitudes of different minimum repulsion energy configura- 
tions. The energetically optimal wavefunction accounts 

for both qualitative ('(-I )' order) and quantitative 

(size of the order parameter) aspects of the supersolid 
phase in the limit of strong nearest neighbour repulsion. 
Our treatment thus represents a remarkable instance in 
which a non-trivial constraint on the lattice scale imposed 
by a dominant term in the Hamiltonian can be taken into 
account in an exact manner in a variational calculation — 
for instance, the celebrated problem of accounting for the 
single-occupancy constraint imposed (by strong coulomb 
repulsion) on holes in the Cu-0 planes of high-Tc super- 
conductors has resisted analytical treatment thus far. 

Our starting point is the observation that the super- 
solid ground state at large V/t must have a wavefunction 
that lies entirely in the subspace of minimum interac- 
tion energy configurations. More precisely, in the classi- 
cal {t — 0) limit, the ground states are diagonal in the 
particle-number basis, and form an extensively degener- 
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FIG. 1: (color online), a) The two possible supersolid states 
at the three-sublattice wavevector Q. Dark (light) links rep- 
resent higher (lower) bond kinetic energy, and sites are color 
coded to reflect mean density, b) Flippable single and double 
hexagons in the dimer representation, and the corresponding 
local density configuration. 

ate set of states corresponding to all minimum repulsion 
energy configurations of particles. To leading order in 
t/V, the slow dynamics induced by the kinetic energy 
term then leads to an effective Hamiltonian Ti. acting in 
within this manifold of states: 

'^^eff^-tJ2vMh+b]b.)T's (2) 

where Vg is the projection operator to the minimum re- 
pulsion energy subspace. 

As H maps to a system of S* = 1/2 spins (5f ~ Ui — 1/2 
where rii is the boson number at site i) with frustrated 
antifcrromagnetic exchange = V between the z com- 
ponents of neighbouring spins, ferromagnetic exchange 
J± = 2t between their x and y components, and magnetic 
field Bz oc fi, the ground states in this t — limit may be 
identified with minimally frustrated states of the classical 
Ising antiferromagnct on the triangular lattice 8]. As is 
well-known, these may be conveniently characterized in 
terms of dimer coverings of the dual honeycomb lattice 
(with a dimer placed on the dual link perpendicular to ev- 
ery frustrated bond of a given spin configuration). In this 
language, Ti-eff is simply a quantum dimer model with 
a ring-exchange term that acts on each pair of adjacent 
hexagons on the dual honeycomb lattice (corresponding 
to each bond {ij) of the triangular lattice): 

^e// = -t J2 (|:«'fe))('3;fe)|+/i.c.) (3) 

Thus, the supersohd behaviour at large V/t should be 
understood in terms of the ground state of this quantum 
dimer model. As was noted in earher literature 0, 
a trial state obtained from an equal amplitude superpo- 
sition of all minimum interaction energy configurations 
(which maps, apart from a global particle-hole transfor- 
mation, to a uniform superposition of dimer covers of the 




FIG. 2: (color online), a) Dimer ensemble with a three- 
sublattice pattern of fugacities which breaks lattice transla- 
tion symmetry at wavevector Q. b) The two types of flippable 
double-hexagons, with two flippable conflgurations each — 
note that light (dark) links correspond to fugacities 1 [z). 

honeycomb lattice) already provides a substantial kinetic 
energy gain, while minimizing the inter-particle repulsion 
by construction. 

Since (5'o|6||'I'o) is clearly proportional to the non- 
zero probability that hexagon i is flippable in the clas- 
sical dimer model, this trial state immediately provides 
a rationale for the persistence of off-diagonal long-range 
order and superfiuidity in the large V/t limit. On the 
other hand, this trial state is unable to fully account for 
the density wave order of the true supersolid state, as 
density correlators in | ^'o) map on to spin correlations of 
the T = classical Ising antiferromagnct on the trian- 
gular lattice which does not support genuine long-range 
order 8], but instead displays power-law order at the 
three-sublattice wavevector Q . 

Here, we focus instead on a variational state obtained 
from a dimer model with two types of links, with different 
dimer fugacities 1 and z > 0, defined on the honeycomb 
lattice in a pattern (Fig[2]) which breaks lattice symmetry 
at wavevector Q for all z 7^ 1: 

l^.ariz)) - ^V^c.W/2(|n+(C,)) + |n_(C,))) (4) 

where Pc^ (z) is the probability of a dimer configuration 
Cd in this dimer model, and the two particle-hole conju- 
gate configurations \n±{Cd)) corresponding to Cd occur 
with the same amplitude in this nodeless wavefunction 
(this follows from the Perron-Frobenius theorem, as was 
noted earlier in a different context [lo|). 

Clearly, |vE'(z)) is characterized by three-sublattice den- 
sity wave order of the (H ) type (see Fig[T]) for z < 1 

while for z > 1, it displays density wave order of the 
(H — 0) type. Furthermore, as the state is constructed 
from a coherent superposition of Fock states with a con- 
siderably wide distribution of average density, |^'(z)) is 
expected to also possess off-diagonal long range order as- 
sociated with superfluidity, at least for z close to z = 1. 

In order to perform an unbiased variational 
study, we need to locate minima of the varia- 
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FIG. 3: (color online). Variational energy per site in units of 
2t as a function of z. 



tional estimate of energy per site (in units of 2t): 
E(z) = {^^ar[z)\'He.ff\^var{z))/2tL'^. To Calculate 
E {z), wc note tha t {^yar{z)\'Pgb\b-jVg\^yar{z)) = 
./PW^~))Pi^^^, where P{^(^^,)) and P(03^,^.p 
are the probabilities that the double hexagon {ij) is in 
the flippahle configurations depicted in their respective 
arguments. This allows us to write 



E{z) - - [2^ P{IC)P{IC) + ^P{2C)P{2C) j (5) 

where P(1C) and P{IC) (P(2C) and P{2C)) are the 
respective probabilities that type 1 (type 2) double- 
hexagons are in flippable configurations C and C (see 
Fig 12). Thus, the variational energy is lowered if the 
expected number of flippable double-hexagons is large. 
Since the number of flippable double hexagons must de- 
crease rapidly to zero in the z ^ as well as the z — > oo 
limit as the dimers freeze into a perfect columnar or pla- 
quette state in these limits, it is immediately clear that 
one or more variational minima of E(z) must occur at 
or in the vicinity of the translationally invariant point 
z = l. 

Calculation of the probabilities P is greatly facilitated 
by the well-known formulation of dimer models on planar 
graphs in terms of Grassmann variables ll|, ll3| : For 
the case at hand, this can be obtained by using the 'arrow 
convention' displayed in Fig[2]to define an antisymmetric 
matrix M, with Mtj = +fJ'{ij) (where the fugacity asso- 
ciated with link (ij)) if an arrow points from point i to j 
and Mij = —^{ij) if the arrow goes from j to i {Mij = 
if i,j are not nearest neighbours, /j.^^^ equals z or 1 as 
shown in Fig [2]). The dimer partition function Z is then 
obtained as Z = \Pf[M]\ = \J[V^I:]exp{S)\, where Pf 
denotes the Pfaffian and S — J2i<j ^hji^ii^j is the ac- 
tion for Grassmann variables tpi defined on sites of the 
honeycomb lattice. 

To proceed further, we represent the honeycomb net 



Q. 

to 



0.12 
0.1 
0.08 
0.06 
0.04 
0.02 




z=0.925 
z=1.075 




-0.025 -0.015 



-0.005 0.005 
5p = p -1/2 



0.015 0.025 



FIG. 4: (color online). Distribution of the total density at 
the global variational minimum at z ~ 0.925, as well as at 
the competing local minimum a.t z ^ 1.07. Data shown is for 
L = 96. 



as an underlying triangular Bravais lattice with Roman 
letter coordinates (representing centers of those hexagons 
whose links all have fugacity z) that is decorated by six 
basis points (corresponding to vertices of such hexagons) 
labeled by Greek letters (Fig [5]), and write the action 
as S" = \Y.x,aY.y,i3^3^,^'^a.x'4^p,y- Transforming to 
Fourier space, we obtain S = ^J2k^a,0 ^^''^'^a,k^f3-k^ 
with the 6x6 matrix M^'^ is given as 

k 



M(fc) = 



R(fc) 
-Rt(fc) 



where is the 3x3 null matrix and R can be written as 



R(/s) 



— z —z e 



z —z 



Next, we note that this Grassmann formulation pro- 
vides a simple prescription for the probability that five 
dimers occupy alternate links /i . . . Z5 on the perimeter of 
double hexagon made up of points z = 1, 2 . . . 10: 



\m=l / 



|(V'lV'2V'3 • •V'lo)| (6) 



where ^.i^ are the fugacities of the 5 links occupied by 
dimers. This allows us to write E{z) as 

E{z) = -2z5|(V'lV'2V'3 ■ -^Z-IO)!! - 2:^|(V'lV'2V'3 • ■V'10)2|(7) 

where the subscripts on the 10-point correlators refer 
to the type of double-hexagon on which they are cal- 
culated (see Fig (2). 



4 



The 10 point correlators involved in this expression 
can be computed in the free field action S using Wick's 
theorem and knowledge of the two point correlators: 



4^ 



expHfc. (a;-y))(-M-i(-fc))„,^(8) 



where the integration is over the BriUouin zone 
e(-7r, tt]. As (ipitpj) = if z and j belong to the 
same sublattice on the honeycomb lattice, one needs to 
keep track of 'only' 5! = 120 contractions in evaluating 
either of the two 10-point correlators. Keeping track of 
these contractions and performing the required k integra- 
tions using MATHEMATICA, we obtain the variational 
energy E{z) shown in Fig[3l As is clear from this figure, 
E(z) has two minima, a global minimum at z ~ 0.9250, 
and local minimum at z ~ 1.0750 with energy only very 
slightly higher (by about w 0.047%) than its value at the 
global minimum. 

Thus, our variational calculation yields a (H ) type 

three sublattice ordered supersolid state in the large V/t 
limit, and suggests that the energy of the (H — 0) su- 
persolid at the same wavevector is only slightly higher 

than that of the (-1 ) supersolid. In order to obtain 

an independent verification of this result, we have also 
used the method of Ref Q to simulate the corresponding 
dimer model and numerically calculate E(z) — the results 
of these calculations are seen to match precisely with the 
results obtained by Grassmann techniques (Fig [3]). 

This simulation also allows us to measure directly the 
solid order parameter -0 — nia + mf,e^'^'/'^ -f- ttIcC^^*/'^, 
where rua, mt, rric are the number densities on the three 
sublattices of the triangular lattice. We obtain \ip\'^ = 
0.03885 for the global minimum dX z — 0.9250, which 
agrees within 10% with the results from QMC [5| ex- 
trapolated to V/t —foo. ( the value in the competing 
(-1- - 0) state is IV'P = 0.03697 at the shghtly higher en- 
ergy local minimum at z = 1.0750). In addition, we also 
measure the histogram of the total density, from which 
one may calculate the ground state expectation values of 
all powers of the density in the supersolid state. While 
the local minimum at slightly higher energy has no spon- 
taneous deviation of density from half-filling, we find that 
the density histogram at the global minimum shows a 
characteristic two-peak structure, reflecting a very small 
(^ 2%) spontaneous deviation from half-filling charac- 
teristic of the (H ) supersolid (Fig[4]). 

To obtain further insight into the smallness of this 
spontaneous density deviation from half-filling, we note 
that reducing z slightly from z = 1 introduces a field 
coupled to the Fourier component of the particle-density 
at wavevector Q. The effect of this field can be under- 
stood by using the well-known 15| height model for- 
mulation of z = 1 dimer model in terms of an action 
Sh = K J (PxiyK)^ for a 'height' field h{x) in terms of 
which one may write the Fourier component of the par- 
ticle density at wavevector Q as ~ exp{inh/3), and 



at wavevector as ptot exp^inh). 

For z < 1 but close to z = 1, the response of the 
density at these two wavevectors can be obtained by cal- 
culating the corresponding susceptibilities, and one finds 
that the response at wavevector Q is much larger than at 
wavevector 0, thereby providing a rationalization for the 
smallness of this symmetry allowed density deviation (in 
comparison with the size of the order parameter) . 

Thus, our variational approach provides a coherent pic- 
ture of the triangular lattice supersolid at large V/t: Ki- 
netic energy effects are seen to select a (-1 ) supersolid 

state with a three-sublattice density wave order parame- 
ter comparable in magnitude to the mean density itself. 
In addition, the competing (H — 0) supersolid is seen to 
be very close in energy to this variational ground state, 
consistent with numerical evidence and analytical argu- 
ments [16i] that the order parameter phase (which distin- 
guishes between these two states) is only weakly pinned 
in the supersolid state. Furthermore, the spontaneous de- 
viation of the density from 1/2 is seen to be a very small 
fraction (~ 2%) of the mean density, consistent with the 
surprisingly small value obtained in earlier QMC studies. 
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